#include <sstream>
#include <iostream>
#include <algorithm>
#include "UPGMA_nearest_neighbor.h"
#include "fill_matrix_lacunae.h"

using std::vector;
using std::string;
using std::cerr;
using std::copy;
using std::ostringstream;
using bcpp_tree::t_tree;
using bcpp_tree::t_node;


namespace upgma_tree
{
//	for (unsigned i = 1; i < cnt_leaves; ++i)
//		for (unsigned j = 0; j < i; ++j)
//	for (unsigned i = 0; i < cnt_leaves - 1 ; ++i)
//		for (unsigned j = i + 1; j < cnt_leaves; ++j)
//	two forms are identical


//----------------------------------------------------------------------------------------

//	upgma_make_tree

//----------------------------------------------------------------------------------------
void upgma_make_tree(	bcpp_tree::t_tree& tree,
						const std::vector<std::string>& names,
						const t_matrix& orig_matrix,
						bool do_fill_lacunae,
						t_progress_indicator& dots)
{
	assert(names.size() == orig_matrix.size());
	unsigned cnt_leaves = names.size();

	// make copy of matrix data so we don't clobber it
	t_matrix matrix(orig_matrix);
	if (do_fill_lacunae)
		fill_lacunae(matrix);



	// prepare a list of clusters each of which starts out with one leaf
	vector<t_node*> cluster;
	for (unsigned i = 0; i < names.size(); ++i)
		cluster.push_back(&tree.add_leaf(names[i], i));


#ifndef NDEBUG
	// make sure matrix is symmetrical!!
	for (unsigned i = 1; i < cnt_leaves; ++i)
		for (unsigned j = 0; j < i; ++j)
			assert(matrix[j][i] == matrix[i][j]);
	for (unsigned i = 0; i < cnt_leaves - 1 ; ++i)
		for (unsigned j = i + 1; j < cnt_leaves; ++j)
			assert(matrix[j][i] == matrix[i][j]);
#endif


	vector<double>	avg			(cnt_leaves, 0.0);
	vector<long> 	cluster_sz	(cnt_leaves, 1);

	// iterate (cluster - 1) times
	unsigned cnt_cluster = cluster.size();
	while (--cnt_cluster)
	{
		++dots;

#ifndef NDEBUG
		// make sure matrix is still symmetrical!!
		for (unsigned i = 0; i < cnt_leaves - 1 ; ++i)
			for (unsigned j = i + 1; j < cnt_leaves; ++j)
				assert(matrix[j][i] == matrix[i][j]);
#endif


		// Find smallest value in matrix
		unsigned	min_i = 0;
		unsigned	min_j = 0;
		double		min_dist = 99999.0;		// ignore dead clusters
		for (unsigned j = 1; j < cnt_leaves; ++j)
		{
			if (!cluster[j])
				continue;
			for (unsigned i = 0; i < j; ++i)
			{
//				if (!cluster[i])
//					continue;
				if (matrix[j][i] < min_dist)
				{
					min_dist = matrix[j][i];
					min_i = i;
					min_j = j;
				}
			}
		}

		double half_dist = matrix[min_j][min_i] * 0.5;

		// merge cluster j into cluster i
		cluster[min_i]->set_length(half_dist - avg[min_i]);
		cluster[min_j]->set_length(half_dist - avg[min_j]);
		cluster[min_i] = &tree.join_nodes(*cluster[min_i], *cluster[min_j]);
		cluster[min_j] = 0;


		// use cluster j distances to re-weight cluster i distances in matrix
		for (unsigned i = 0; i < cnt_leaves; ++i)
		{
			if (cluster[i] &&
				matrix[i][min_i] + matrix[i][min_j] < 5000.0)
			{
				matrix[i][min_i] =
				matrix[min_i][i] =
									(matrix[i][min_i] * cluster_sz[min_i] +
									 matrix[i][min_j] * cluster_sz[min_j]	)
									 /
									 (cluster_sz[min_i] + cluster_sz[min_j]);
			}
		}

		// cluster j distances are discarded from matrix
		for (unsigned i = 0; i < cnt_leaves; ++i)
			matrix[i][min_j] = 100000.0;
		for (unsigned i = 0; i < cnt_leaves; ++i)
			matrix[min_j][i] = 100000.0;

		// update cluster i size with cluster j
		cluster_sz[min_i] += cluster_sz[min_j];
		cluster_sz[min_j] = 0;
		avg[min_i] = half_dist;
	}

	// last standing cluster still attached to root
	tree.finish();
	return;
}

//________________________________________________________________________________________

//		upgma_make_tree

//			wrapper

//________________________________________________________________________________________
void upgma_make_tree(	bcpp_tree::t_tree& tree,
						const std::vector<std::string>& names,
						t_matrix& matrix,
						bool fill_lacunae)
{
	ostringstream ignore_strm;
	t_progress_indicator dots(ignore_strm, 10000);
	upgma_make_tree(tree, names, matrix, fill_lacunae, dots);
}




}

